PDE DSA

This notebook computes the numerical simulation of the PDE system
$\left(\partial_t+\partial_s\right)y_S = -\left(\lambda+\int^\infty_0\beta y_Idu\right)y_S$,
$\left(\partial_t+\partial_s\right)y_V = -\left(\alpha+(1-\alpha)\int^\infty_0\beta y_Idu\right)y_V$,
$\left(\partial_t+\partial_s\right)y_V = -\gamma y_I$.
Above the equation coefficients are functions of (t,s) as described in the main paper. Additionally, the system has the boundary data
$y_S(0,s) = f_S(s)$ and $y_S(t,0)=0$ for all $t\geq 0$,
$y_V(0,s) = 0$ for all $s\geq 0$ and $y_V(t,0) = \int^\infty_0\lambda y_S du$,
$y_I(0,s) = \rho f_I(s)$ and $y_I(t,0) = \int^\infty_0y_S\int^\infty_0\beta y_Iduds + \int^\infty_0\left(1-\alpha\right)y_V\int^\infty_0\beta y_Iduds$.

The user specifies the parameter values for these equation coefficients in the data() function of the parameters.jl file. Then one runs this notebook.
NOTE: One does not edit the data!() function.

include("auxilliary.jl"); include("parameters.jl"); include("flow.jl"); include("abcmc.jl");

Simulate the system

prm = data();
ysol,yʳsol=pdesolve(;prm=prm);

Simulation progress: 0.015625/1
Simulation progress: 0.03125/1
Simulation progress: 0.046875/1
Simulation progress: 0.0625/1
Simulation progress: 0.078125/1
Simulation progress: 0.09375/1
Simulation progress: 0.109375/1
Simulation progress: 0.125/1
Simulation progress: 0.140625/1
Simulation progress: 0.15625/1
Simulation progress: 0.171875/1
Simulation progress: 0.1875/1
Simulation progress: 0.203125/1
Simulation progress: 0.21875/1
Simulation progress: 0.234375/1
Simulation progress: 0.25/1
Simulation progress: 0.265625/1
Simulation progress: 0.28125/1
Simulation progress: 0.296875/1
Simulation progress: 0.3125/1
Simulation progress: 0.328125/1
Simulation progress: 0.34375/1
Simulation progress: 0.359375/1
Simulation progress: 0.375/1
Simulation progress: 0.390625/1
Simulation progress: 0.40625/1
Simulation progress: 0.421875/1
Simulation progress: 0.4375/1
Simulation progress: 0.453125/1
Simulation progress: 0.46875/1
Simulation progress: 0.484375/1
Simulation progress: 0.5/1
Simulation progress: 0.515625/1
Simulation progress: 0.53125/1
Simulation progress: 0.546875/1
Simulation progress: 0.5625/1
Simulation progress: 0.578125/1
Simulation progress: 0.59375/1
Simulation progress: 0.609375/1
Simulation progress: 0.625/1
Simulation progress: 0.640625/1
Simulation progress: 0.65625/1
Simulation progress: 0.671875/1
Simulation progress: 0.6875/1
Simulation progress: 0.703125/1
Simulation progress: 0.71875/1
Simulation progress: 0.734375/1
Simulation progress: 0.75/1
Simulation progress: 0.765625/1
Simulation progress: 0.78125/1
Simulation progress: 0.796875/1
Simulation progress: 0.8125/1
Simulation progress: 0.828125/1
Simulation progress: 0.84375/1
Simulation progress: 0.859375/1
Simulation progress: 0.875/1
Simulation progress: 0.890625/1
Simulation progress: 0.90625/1
Simulation progress: 0.921875/1
Simulation progress: 0.9375/1
Simulation progress: 0.953125/1
Simulation progress: 0.96875/1
Simulation progress: 0.984375/1
Simulation progress: 1.0/1
Number of times Δt was too small for refinement: 0

Plot equation coefficients

Below we plot the equation coefficients used in the simulation according to the parameters.jl file. Legends labeled true or discrete refers respectively to the curve given by its analytic expression versus the approximation obtained by its nodal interpolant spline. One expects these approximations to coincide for suitably fine meshes. The abbreviation cap is only significant when one the two Weibull distributions for either the contact interval $\beta$ and the recovery period $\gamma$ have largely disjoint supports with wide separation. Because then one of the hazards is exceedingly large over the domain of simulation, this can lead to numerical instability so it can be convenient to cap the hazard to a maximum value. In the default simulation parameters, one sees that this was not used as the curves overlap.
NOTE: The symbol $\partial$ in legend titles means the equation coefficient's $\left(\partial_t+\partial_s\right)$ derivative.

plot(:α;prm=prm);plot!(size=(600,300))

plot(:β;prm=prm);plot!(size=(600,300))

plot(:γ;prm=prm);plot!(size=(600,300))

plot(:Weibull;prm=prm);plot!(size=(450,300),titlefontsize=12)

plot(:fˢ;prm=prm);plot!(size=(450,300))

plot(:fⁱ;prm=prm);plot!(size=(450,300))

Plot the solution

plot(ysol);plot!(size=(900,300))

plot(ysol,yʳsol;prm=prm);plot!(size=(450,300))